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ABSTRACT 

Both simulations and observations indicate that stars form in filamentary, hierarchically clustered 
associations, most of which disperse into their galactic field once feedback destroys their parent clouds. 
However, during their early evolution in these substructured environments, stars can undergo close 
encounters with one another that might have significant impacts on their protoplanetary disks or young 
planetary systems. We perform N-body simulations of the early evolution of dissolving, substructured 
clusters with a wide range of properties, with the aim of quantifying the expected number and orbital 
element distributions of encounters as a function of cluster properties. We show that the presence 
of substructure both boosts the encounter rate and modifies the distribution of encounter velocities 
compared to what would be expected for a dynamically relaxed cluster. However, the boost only lasts 
for a dynamical time, and as a result the overall number of encounters expected remains low enough 
that gravitational stripping is unlikely to be a significant effect for the vast majority of star-forming 
environments in the Galaxy. We briefly discuss the implications of this result for models of the origin 
of the Solar System, and of free-floating planets. We also provide tabulated encounter rates and 
orbital element distributions suitable for inclusion in population synthesis models of planet formation 
in a clustered environment. 

Subject headings: open clusters and associations: general — planets and satellites: formation — stars: 
kinematics and dynamics 



1. INTRODUCTION 



Stars 



that 
tend 



form in giant molecular clouds (GMCs) 
possess a high degree of subst ructure. They tena 
to b e clumpy and filamentary (jWilliams et al.l 11994 
120001 ) ■ almost certai nly as a res i ilt of pervasive su- 
perso nic turbulence ()Larsonlll981l : iMac Low fc KlessenI 
l2004f ). Stars that form out of these clouds in- 
herit the substructures of the ir parents, le ading 
to a hierarchy of clustering (iLada fc Ladal 120031 : 
iBressert et"alll2010t ICutermuth et al.ll201lD , and to self- 
simila r (fractal) s tructures within star clust ers (iLarsonl 
19951: lElmegreenI l2000t lElm egreen~fc Elmcgreed |2001t 
Cartwright fc WhitworthI 12004: Chen et al. 2005). Nu- 
merical simulations o f star formation produce simi- 
lar results (Klcssen fc Burkertl[20 00: Bon nell et al.ll2003l : 
lOffner ct al. 2009; Krui nholz et al.ll2012l) . 

Aarseth fc Hills (1974) were the first to study the evo- 
lution of star clusters with initial substructure. They 
found that any substructure initially present was typi- 
cally destroyed within one free- fall time, and observations 
gener ally supp ort this pictu re (jCartwright fc WhitworthI 
120041 : Schmeia et al.l 120081 ). However, simulations in- 
dicate that substructure can be destroyed quickly 
only in initially subvirial clusters, a case character- 
ized by an initial collapse of the star cluster towards 
the center of mass followed by a chaotic evolution- 
ary p hase ([Goodwin fc WhitworthI 12004 I Allison et"all 
|2010( ). For supervirial clusters, on the other hand, sub- 
structure can survive for up to several crossing times 
(jGoodwin fc Whitworthll2004 ). 

Both observations and theory suggest that clusters typ- 
ically form subvirially with respect to the gas, though not 
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necessarily with respect to the stars alone (jFtiresz et al.l 
2008; Tobin et al. 2009; Offner ct al. 2009). However, the 
star formation process is inefficient, with relatively small 
amounts of the mass in a given molecular cloud being 
converted into stars, which then expel the remainder of 
the cloud into back into the diffuse interstellar medium 
throu gh their rad i ation, winds, and supernovac (Hilla 
1980; Lada 1999; Lada fc Lada 2003; Matzner 2002 
Krumholz et al 2006 : FaU et al. 2010: Goldbaum et atl 
2011; Krui issenll2012D . Once the gas is expelled, stars 
disperse into the field, with only a minority remaining in 
bound clusters for many dynamical times after gas dis- 
persal. As a result, even if stars are born subvirial with 
respect to the gas, they may be rendered supervirial by 
its rapid dispersal. Real star clusters may therefore ex- 
perience periods of both subvirial and supervirial evolu- 
tionQ 

Whether subvirial or supervirial, this early evolution- 
ary stage is of considerable interest for the problem of 
planet formation. In denser environments and mas- 
sive clusters containing massive sta rs, where a signifi- 
cant fraction of stars appear to form (|Lada fc Ladall200"3l : 
iChandar et al.ll20Tol ). close passages between Solar- type 
stars and massive stars may lead to the photoevaporation 
of protoplanetary disk s and mod ification of the planet 
form a tion process (Ada ms et al.l 2004 ; Throop fc Ballvl 
120051 : lAdamd [20101 and references therein). Close en- 
counters with passing stars can also gravitationally dis- 

^ A brief comment on terminology: some authors use the word 
"cluster" to refer to any significant stellar over density regardless of 
its dynamical state, while others use the term to refer exclusively 
to stellar structures that remain bound after gas dispersal. We 
follow the former approach, and refer to our objects as clusters 
even though they are unbound. 



rupt both disks and planetary systems, potentially trun- 
cating disks, exciting planetary orbits, or ejecting plan- 
ets completely. Such encounters have been suggested as a 
potential explanation for such diverse o bservations as the 
existence o f free-floating planets (e.g. iSumi et al.l 120111 : 
iVeras fc: Raymond 2012) and the s tructure of the Kuiper 
Belt (e.g. iLestrade et al.l 120 lit iJimenez-Torres et al.l 
120111 ). The need to avoid disruptive encounters has also 
been used as a co nstraint on the potential birth environ- 
ment of the Sun (lAdams fc Laughlinl l200ll : lAdams et al.l 
I2OO6) . Our Solar System is remarkably well ordered com- 
pared to many extrasolar planetary systems, with all of 
the planets on nearly circular orbits (every planet ex- 
cept Mercury has e < 0.09) , while the Kuiper belt is 
also relative l y undisturbed. lAdams fc Laughliiil (|2001[ ) 
and lAdamsl J2010) have argued that this implies that 
the Sun could have been formed in a cluster no larger 
than ^ lO'^ stars, though this conclus i on ha s recently 
been questioned bv lDukes fc Krumhold ()2012[ ). 

A crucial input to all these questions is the rate and 
distribution of orbital elements of the encounters that 
a star in a cluster will experience. These are a neces- 
sary ingredient for population synthesis mo dels for planet 
formation in clustered environments (e .g. | Adamsl 120101 : 
IDukes fc Krumholzl [Ml lOvelar et al.l |2012|). While a 
number of auth ors have measured these distributions 
numerically (e. g. 'Bonnell ct al."2001'; 'Adams et alll2006l: 
ISpurzem etaLi 2009; Olczak ct al. 2006, 2010), none thus 
far have done so in the context of a dispersing, initially 
highly-substructured cluster, which modern observations 
suggest is the typical condition for the formation of most 
stars. 

Several authors have recently studied properties of ini- 
tially substructurcd (fractal) clusters in slightly different 
contexts. Parker ct al. (2011) find that reproducing the 
observed present-day binary properties of the ONC re- 
quire that it have formed with a high de gree of substruc- 
ture a nd a high initial binary fraction. iParker fc Meveii 
(|2012D find that the surface density of fractal star clusters 
decreases rapidly in time, which implies that a large frac- 
tion of star- star encounte r s will occur very early on in the 
cluster life. I Smith et al.l (|2013l ) found that gas removal 
after multiple crossing times typically results in relaxed 
clusters, with no remaining substructure, whereas quick 
gas expulsion before one crossing time leads to a stochas- 
tic, unpredictable outcome. 

T he two papers closest to our work are lAdams et al.l 
()2006J) . who calculate encounter rates in both dispersing 
and cold clust ers, but do not include any initial sub- 
structure, and IParker fc Quanzl (|2012D who do include 
substructure and study how star clusters affects orbital 
elements of planetary systems. We add to these studies 
by conducting a series of N-body simulations of dispers- 
ing, fractal star clusters across a much broader param- 
eter space than has been considered before. We con- 
sider a wide range of dynamical environments, from un- 
bound, supervirial stellar associations to subvirial stars 
in a gas-dominated clump that are subsequently unbound 
by gas expulsion. For each simulation we track every 
event in which two or three stars pass within 1000 AU 
of one another, giving a nearly complete dynamic pro- 
file o£_£OSsibleinteractions. Our work expan d s on that 
of lAdams et all ()2006t ) and IParker fc Quanzl (|2012D by 
surveying a significantly broader parameter space, with 



cluster masses from 30 — 30, 000 Mq and surface den- 
sities from 0.1 — 3.0 g cm~^. This broad survey allows 
us to measure how the results depend on cluster mass 
and surface density, and thereby to extrapolate into the 
regime of high mass and surface density clusters that are 
too computationally expensive to simulate directly. 

The remainder of this paper is as follows. Section [5] 
discusses the model parameters, the initial conditions 
for the clusters, and the simulations and data reduction 
methods, and defines the statistical distributions of inter- 
est. Section [3] details our results, and Section [4] discusses 
their implications. Our conclusions are presented in sec- 
tion [S] 

2. METHODS 

To study stellar encounters in dissolving clusters, we 
perform an ensemble of N-body simulations using a 
modified version of the numerical integrator NB0DY6 
(jAarsethI [1999) . Below we describe the parameters in 
our simulations, the initial conditions, and how we pro- 
cess the resulting data. 

2.1. Simulation Parameters and Initial Conditions 

We characterize clusters by four parameters: the virial 
ratio, Q, defined as the ratio of kinetic to potential en- 
ergy (so that a cluster in equilibrium has Q = 0.5), the 
fractal dimension D, the stellar mass Mc, and the clus- 
ter surface density, S^. We describe below how we use 
these parameters to set up the initial conditions in our 
simulations. We consider four combinations of Q and D, 
and for each combination we then simulate clusters with 
a broad range of masses Mc and surface densities Sc. 

We use the surface density Sc rather than the ra- 
dius Re or the volume density pc as a parameter for 
two reasons. First, while the volume density determines 
encounter rates, the quantity of interest for the stand- 
point of studying how clusters affect planetary systems 
is the total number of encounters a star can expect to 
experience over the cluster lifetime, not the encounter 
rate. The natural time scale for a disrupting cluster is 
the crossing time, and the total number of encounters 
per crossing time depen ds on the surface density rather 
than the volume density ([Dukes fc KrumhoI3l2012l ). Sec- 
ond, observations of cluster-forming gas clumps appear 
to indicate that, while clusters span a very wide range of 
volume densities and radii, the y form a sequence of rela- 
tively constant surface density (I Fall e t al. 2010,, and ref- 
erences therein). Thus in discussing embedded clusters 
that have until recently been dominated by the potential 
of the gas, it is also natural to work in terms of surface 
rather than volume density. 

Our base case is a cluster with Q — 0.75 and D — 2.2; 
this value of virial ratio corresponds approximately 
to a cluster that has just expelled its residual gas 
but has not been completely unbound by the process, 
and the fractal dimension describes a cluster with 
a moderate degree of substructure. This is consis- 
tent with observations which have typically found 
that D goe s from 1.9 to 2.5 (Falgarone ct al. 199J; 
Vogelaar fc Wakker 1994; Elmcgrcen & Falgarone 199^; 
de La Fuente Marcos fc de La Fucnte Mar cos ,2006 : 
Sanchez et al.l |2010[ ). These runs generally result in 



majority of the stars escaping promptly, but some 
remaining as a bound structure for long times. Our 



second case is Q — 1.25 and D — 2.2, corresponding 
to a cluster that has been completely unbound by gas 
expulsion, or that was never bound in the first place. 
In these runs essentially all the stars disperse. Our 
third case is a model with D = 1.6 with Q = 0.75, 
corresponding to a case like the base model but with 
more substructure. In our fourth model, we explicitly 
include a phase in which the stars are confined by an 
external potential, which we rapidly remove after four 
crossing times, where the crossing time is 
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(1) 



Our motivation for this choice is that observations indi- 
cate that the lifetime of the embedded phase o f star clus- 
ter formation is roughly four crossing times (jTan et al.l 
1200 (f). or possibly even less (Elmcgrcen 2000l). We should 
note here that substructure is typically eras ed after sev- 
eral c rossing times in a confined potential (jSmith et al.l 
l2013f ). The first three cases assume that the gas is ex- 
pelled very early while the substructure still remains. 
While it is present, we describe the gas potential by a 
Plummer model, 
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and we choose Mgas = (0.7/0.3)Mc, so that the gas mass 
is 70% of the total gas plus stellar mass. This cluster has 
Q = 0.3 and D = 2.2 initially. We should note that this 
value of the virial ratio is computed using only the po- 
tential energy due to the interactions between stars, not 
the coupling of the stars to the gas. The total potential 
energy is 
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where N is the number of stars. 



'■I 1=1 

(3) 
is the distance be- 
tween the i*'' and j"^ stars, and r^ is the radial position 
of the i^^ star. This means that the real virial ratio is 
then T jU tot- The gas mass term dominates, which leaves 
us with Q < 0.1, an extremely subvirial case. This effec- 
tively gives us a bound on how important the effect of 
gas might be on the evolution of the cluster. 

We summarize the model parameters in Table [I] and 
the number of independent realizations we perform at 
each {Mc, Sc) combination in Tables |3] - El The num- 
bers of runs for each {Ale, ^c) value are chosen so that 
the number of interactions, and thus the statistical error 
on our results, is roughly constant. This implies a large 
number of runs for small, low surface density cases, and 
a smaller number of runs for more massive, higher sur- 
face density cases. As we will see below when we discuss 
our error budget, this does limit our accuracy to some 
extent in the high mass and surface density regime. Un- 
fortunately this regime is too computationally-costly to 
allow a significantly larger number of simulations. The 
relatively small number of simulations limits our accu- 
racy, but even with this limitation we show below that 
the errors on our measured encounter rates are typically 
no more than ~ 10%. 



TABLE 1 
Model parameters 



Name 


Q 


D Gas? 


Q0.75D2.2 


0.75 


2.2 No 


Q1.25D2.2 


1.25 


2.2 No 


Q0.75D1.6 


0.75 


1.6 No 


Gas 


0.3 


2.2 Yes 


TABLE 2 




Cluster mass and number 


OF 


STARS, 




log(M,/M0) N 


= Mc/m 


1.5 




54 


2.0 




170 


2.5 




540 


3.0 




1707 


3.5 




5400 


4.0 




17077 


4.5 




54003 



TABLE 3 
Number of realizations for model Q0.75D2.2 



Sc 






logiMc/Me) 






[g cm-2] 


1.5 


2.0 


2.5 


3.0 


3.5 


4.0 


4.5 


3.0 


600 


150 


40 


8 


5 


2 





1.0 


1000 


200 


50 


15 


5 


4 


3 


0.5 


1500 


300 


75 


20 


10 


5 


3 


0.1 


4000 


800 


150 


40 


20 


5 


4 



TABLE 4 

Number of realizations for model 

Q1.25D2.2 



Se 




log(M,/M0) 




[g cm-2] 


1.5 


2.0 


2.5 


3.0 


3.5 


3.0 


300 


100 


20 


10 


3 


1.0 


500 


175 


30 


15 


4 


0.5 


700 


175 


40 


15 


4 


0.1 


1000 


275 


75 


30 


10 



TABLE 5 

Number of realizations for model 
Q0.75D1.6 



[g cm-2] 


1.5 


log(M,/M0) 
2.0 2.5 3.0 


3.5 


1.0 
0.5 
0.1 


500 
500 
500 


100 15 5 
100 30 6 
150 40 8 




5 



TABLE 6 
Number of realizations for model Gas 



Sc 




log(M,/M0) 






[g cm-2] 


1.5 


2.0 


2.5 


3.0 


3.5 


4.0 


3.0 


175 


50 


15 


6 


3 





1.0 


225 


75 


25 


8 


4 


2 


0.5 


275 


100 


30 


10 


5 


3 


0.1 


350 


125 


40 


12 


6 


4 



2.2. Initial Conditions 

We initialize our clusters using the fr actal initial con- 
ditions model with sliRht modifications (jScallv fc Clarkd 
l2002tlGoodwin fc Whitworth.2061 . We refer the reader 
to the second paper for full details of the method, which 
we briefly summarize below. To generate the cluster we 
start by defining a cube with sides of length 2 (in arbi- 
trary units, which will be scaled later to give the cor- 
rect physical units), centered at the origin. This cube is 
subdivided into the 8 Cartesian sectors, and a first gen- 
eration particle is placed at the center of each subcube. 
Each of these first generation cubes is then subdivided 
again, with second generation particles placed at the cen- 
ter of each second generation subcube. We repeat this 
subdivision procedure, with one additional constraint for 
the second and subsequent generations: each parent par- 
ticle only has a probability 2^~^ of producing offspring. 
When D = 3 this ensures that all positions are equally 
populated and there is no substructure, but for D < S 
parts of the cube will be empty, yielding substructure. At 
each generation g, we also add a random displacement of 
position of magnitude 2~^~^ to prevent the development 
of an overly gridded structure. 

We repeat this procedure until the number of particles 
generated greatly exceeds the number we will actually 
use in the simulation. We then randomly select a subset 
of the points with radius less than 1 (in our arbitrary 
units) to be the initial locations of our stars. The radial 
positions of the stars are then multiplied by a factor 



Me 
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(4) 



so that the average surface density of the cluster is Ec. 

The number of stars is simply Mc/rh, where fh is the 
mean stellar mass for our chosen IMF (see below). Ta- 
ble [5] gives the correspondence between the cluster mass 
and the number of stars N. Note that this means that 
for a given Mc the actual cluster mass may be slightly 
larger or smaller, depending on drawing from the IMF. 
We therefore interpret Mc as the expectation value of the 
cluster mass, though deviations from this value are small 
as long as Mc ^ to. 

We assign initial velocities to the stars using a recur- 
sive procedure to ensure that positions and velocities are 
correlated, as suggested by observations and simulations. 
At each generation we assign a random scalar velocity 
drawn from a Maxwellian distribution to each particle 



p{v,g) oc w^ exp —2^^ 



2<o 



(5) 



The direction of the velocity vector is chosen randomly. 
A particle's velocity is the value produced by this draw- 
ing added to the velocity of its parent. Since the mag- 
nitude of the velocity perturbation decays with genera- 
tion, the positions of the stars are then highly dependent 
on the velocities of the first few parents. Note that the 
choice of a^ g is arbitrary, since we scale the final speeds 
so that the cluster has a specified virial ratio (see below). 
Figured] shows an example of a cluster generated via this 
procedure. 
Finally, we assign stellar masses by randomly drawing 




-0.5 0,0 0.5 

X [N-body units] 



Fig. 1. — Asterisks indicate the positions of stars in an example 
cluster projected onto the xy-plane; colors indicate stars' z veloci- 
ties. Notice that velocities and positions are correlated. 

from an extended version of the IKroupal (j2002[) IMF, 



p{m) ex ^ m^^'^, 

,-2.3 



0.08 < m/Mg < 0.1 

0.1 < tu/Mq < 0.5 . (6) 

0.5 < to/Mq < 120 



This IMF yields a mean stellar mass to sa 0.59 Mq. Once 
we have drawn all the stellar masses, positions, and ve- 
locities, we scale the velocities by a constant factor so 
that the initial virial ratio is the desired value. 

2.3. Simulations and Analysis 

We simulate the evolution of each cluster for 5 cross- 
ing times, except for the Gas runs, which we compute 
for 9 crossing times (i.e. 5 crossing times after the gas 
potential is removed). Since we are interested in close 
encounters for Solar-like stars, we track every instance 
in which a star of mass 0.8 — 1.2 Mq passes within 1000 
AU of another star and the pair has a non-negative cen- 
ter of mass energy; the latter condition excludes cases 
where the two stars form a binary. In addition, once we 
have found a pair of stars that meet this criterion, we 
also check for any other star passing within 1000 AU of 
either body. 

The raw data we obtain from NB0DY6 is a list of 
2-body and 3-body interactions with positions, masses, 
velocities, indices and the time. Almost all interactions 
are recorded as a time series, since the stars involved are 
within 1000 AU of one another for more than a single 
simulation time step. Our end goal is to use these time 
series to calculate the distribution of impact parameters 
b and relative velocities at infinity Uoo for encounters in 
clusters. For a given pair of particle positions and ve- 
locities, it is straightforward to compute what values of 
b and Woo would be required to produce that particular 
separation and relative velocity. However, in practice 
interactions are often complex, particularly in the high 
surface density cases where stars are tightly packed, and 
the values of b and Woo that one computes in this manner 
are not constant over the time series of positions and ve- 
locities that describes a particular interaction. For this 
reason, we must first classify interactions in order to de- 
cide how to analyze them. Our classification scheme is 
as follows. 

1+1 interactions — These are true single star-by-single 
star scattering events. Two bodies are said to be well- 
described as a 1-1-1 interaction if, for that pair of parti- 



cles, the set of impact parameters bk and relative veloc- 
ities Voc,k we conipute from our time series satisfy the 
condition that 



CTb 



<Ti_ 



< Ti+i, 



(7) 



where ab and b are the standard deviation and mean of 
the time series bk and similarly for CTi,^ and Woo, and 
Ti_|_i = 0.1 is the tolerance ratio we adopt for 1+1 inter- 
actions. This value is somewhat arbitrary, but provides 
a reasonable separation between cases where two inter- 
acting stars have their orbits perturbed slightly by the 
potential of other stars during the interaction, and cases 
where another nearby star provides a large perturbation 
to the orbits. For 1-1-1 interactions, we record b and Voo 
as the impact parameter and relative velocity of the en- 
counter. 

1+2 interactions — These are events in which a single 
star scatters off a binary. We label an encounter involv- 
ing three stars within 1000 AU of one another as a 1-1-2 
interaction if two conditions are met. First, exactly one 
pair of the three stars must be gravitationally bound, 
while the other pairs are unbound. Second, if we re- 
place this gravitationally bound pair by a single star at 
the center of mass position and velocity of the pair, and 
with a mass equal to the sum of the pair's masses, and 
we compute a set of impact parameters bk and relative 
velocities Voo.k between this binary and the remaining 
star, we find that 



o-fe 



<Ti_ 






<T, 



1+2, 



(8) 



where T1+2 — 0.3. We set the tolerance somewhat higher 
than in the 1-1-1 case because, even in the absence of 
perturbations from external stars, tidal forces exerted by 
the binary on the single star may lead to some exchange 
of energy and angular momentum between the binary's 
internal energy and angular momentum and that of the 
orbit of the binary and the single _star about one another. 
For 1-1-2 interactions, we record b and Voo as the impact 
parameter and relative velocity of the encounter. 

1+1+1 interactions — These are events in which three 
unbound stars encounter one another. We classify an 
event as 1-1-1 -1-1 if no pair composed of two of the three 
stars involved is mutually gravitationally bound. We de- 
compose encounters of this type into three 1-1-1 events; 
since these 1-1-1 events clearly will not satisfy the toler- 
ance criteria for 1-1-1 events, we simply calculate b and 
Uoo in these cases using the positions and velocities the 
stellar pairs have at their point of closest approach. 

Complex interactions — These are events which do not fall 
into one of the above categories. They may, for example, 
be cases where a metastable hierarchical multiple star 
systems forms and then dissolves some time later. These 
interactions do not have well-defined orbital elements. 
We do not attempt to define an impact parameter or 
relative velocity in these cases, and we do not include 
them in our statistical distributions of b and Uqo. We 
do, however, record such interactions and include them 
in our total counts of events. 



2.4. Statistical Distributions 

For each set of simulations, we are interested in three 
quantities. The first is simply the expected number of 
encounters A'onc within 1000 AU. The other two are the 
distributions of impact parameters p{b) and relative ve- 
locities p{voo) that describe these encounters. For a fully 
relaxed cluster, we expect these to follow 



p{b) ex b, 



and 



p{v^)(xvi^exp[ --^ 



(9) 
(10) 



but they need not for fractal, dispersing clusters that 
have not had time to relax. To evaluate the distribu- 
tions from our simulations, we bin all our encounters into 
A'bin = 20 equally-spaced bins of impact parameter from 
0— 1000 AU, and into A'bin equally-spaced bins of relative 
velocity from — 20 km s^^. 

3. RESULTS 

3.1. Base Case (Q0.75D2.2) 

We first describe the results of our base case, model 
Q0.75D2.2, and then in subsequent sections describe how 
the results change for the other models. We report the 
quantitative results for all models in Tables [7] - [TOl 

3.1.1. Distributions of Encounter Velocities and Impact 
Parameters 

Figure [2] shows the distributions of impact parameter 
for two example runs, one at low and one at high surface 
density. We find that the distribution of impact parame- 
ters follows equation Q, p{b) oc &, very closely for many 
of our cases, and in all cases the distribution of 1 -I- 1 
events follows a linear trend. We find a deviation from 
linearity with the 1 -f 1 + 1 events. This is not terribly 
surprising since we have made a rather large assumption 
that we can reliably describe a 3-body event as three 2- 
body events. As 3-body interactions become more preva- 
lent with higher mass (for reasons to be discussed later) 
the deviation of the overall distribution from linearity 
tends to increase (blue in the figures). However, even for 
the highest surface density cases we consider, the overall 
distribution of impact parameters summed over all 1-1-1, 
1 + 2, and 1 + 1 + 1 events remains reasonably linear. 

While the distribution of impact parameters is close 
to what one would expect for a relaxed cluster, we find 
that the distribution of relative velocities is strongly non- 
Maxwellian in all our simulations. We show some exam- 
ples of the distributions p{voo) from our simulations in 
Figure [31 The deviation from Maxwellian is not sur- 
prising, given the correlated position-velocity distribu- 
tion with which we begin, and that is observed in young 
clusters. 

3.1.2. Number of Encounters 

We now turn to the question of the number of encoun- 
ters and their typical velocities as a function of Mc and 
Ec. First note that nearly all of the events for these clus- 
ters occur in the first crossing time. Figure U shows an 
example of the temporal distribution of encounters in one 
of our cases; other combinations of mass and surface den- 
sity are similar or even more heavily weighted toward en- 
counters occurring during the first crossing time. These 



TABLE 7 
Encounter Statistics for model Q0.75D2.2 



log(Mc/Af0) 


[g cm-2] 


A''cnc 


^Poisson 


•^sample 


^.median 

[km s-i] 


Ur 


1.5 


0.1 


0.85 


0.01 


0.00 


1.88 


0.01 


1.5 


0.5 


3.82 


0.03 


0.00 


2.71 


0.01 


1.5 


1.0 


7.08 


0.04 


0.01 


3.20 


0.01 


1.5 


3.0 


9.64 


0.07 


0.01 


4.23 


0.01 


2.0 


0.1 


1.78 


0.01 


0.00 


2.48 


0.01 


2.0 


0.5 


8.25 


0.05 


0.02 


3.62 


0.01 


2.0 


1.0 


13.69 


0.08 


0.04 


4.18 


0.01 


2.0 


3.0 


30.21 


0.14 


0.10 


5.25 


0.01 


2.5 


0.1 


2.87 


0.02 


0.01 


3.06 


0.01 


2.5 


0.5 


13.63 


0.07 


0.08 


4.60 


0.01 


2.5 


1.0 


23.83 


0.12 


0.20 


5.13 


0.01 


2.5 


3.0 


76.39 


0.23 


0.89 


6.79 


0.01 


3.0 


0.1 


3.46 


0.03 


0.04 


3.41 


0.01 


3.0 


0.5 


17.92 


0.09 


0.54 


5.34 


0.03 


3.0 


1.0 


34.59 


0.15 


0.77 


6.25 


0.02 


3.0 


3.0 


129.41 


0.38 


10.25 


7.37 


0.08 


3.5 


0.1 


4.82 


0.03 


0.06 


4.02 


0.01 


3.5 


0.5 


31.00 


0.10 


1.22 


6.65 


0.04 


3.5 


1.0 


77.44 


0.22 


4.58 


7.58 


0.06 


3.5 


3.0 


160.27 


0.31 


6.49 


9.92 


0.04 


4.0 


0.1 


5.72 


0.03 


0.42 


4.35 


0.07 


4.0 


0.5 


47.44 


0.09 


3.88 


7.72 


0.08 


4.0 


1.0 


92.48 


0.15 


8.95 


9.09 


0.10 


4.0 


3.0 


367.63 


0.43 


60.44 


11.59 


0.16 


4.5 


0.1 


9.01 


0.03 


0.10 


5.47 


0.01 


4.5 


0.5 


46.39 


0.07 


3.65 


8.29 


0.08 


4.5 


1.0 


71.23 


0.09 


10.13 


8.87 


0.14 



Note. — A^cnc is the mean number of encounters within 1000 AU per 
star for a Sun-Ukc star over the full duration of the simulation, o-ptjisson 
and CTgamplc arc the Poisson and parameter space sampling errors on JVenc, 
and (Tr is the total relative error <5 A^onc /A^enc considering both sources; see 
Equations I ll7t . II18I I. and I ll9t . t)m°dian jg ^j^g median encounter velocity. 



TABLE 8 
Encounter Statistics for model Q1.25D2.2 



log(Mc/M0) 


[g cm-2] 


A^enc 


^Poisson 


(^saniplc 


„, median 

[km s-i] 


CTr 


1.5 


0.1 


0.47 


0.01 


0.00 


1.85 


0.02 


1.5 


0.5 


2.04 


0.03 


0.00 


2.73 


0.01 


1.5 


1.0 


3.55 


0.04 


0.01 


3.17 


0.01 


1.5 


3.0 


9.93 


0.10 


0.03 


4.30 


0.01 


2.0 


0.1 


1.08 


0.02 


0.00 


2.31 


0.02 


2.0 


0.5 


5.34 


0.05 


0.02 


3.51 


0.01 


2.0 


1.0 


8.50 


0.07 


0.03 


4.17 


0.01 


2.0 


3.0 


21.40 


0.13 


0.10 


5.19 


0.01 


2.5 


0.1 


2.27 


0.03 


0.02 


2.93 


0.02 


2.5 


0.5 


8.31 


0.08 


0.09 


3.65 


0.01 


2.5 


1.0 


14.05 


0.12 


0.38 


4.80 


0.03 


2.5 


3.0 


63.96 


0.30 


1.44 


6.52 


0.02 


3.0 


0.1 


2.16 


0.03 


0.03 


3.10 


0.02 


3.0 


0.5 


9.31 


0.08 


0.19 


4.65 


0.02 


3.0 


1.0 


28.44 


0.13 


1.35 


5.53 


0.05 


3.0 


3.0 


68.52 


0.25 


3.64 


7.16 


0.05 


3.5 


0.1 


3.30 


0.03 


0.08 


3.70 


0.03 


3.5 


0.5 


22.07 


0.13 


2.62 


5.79 


0.12 


3.5 


1.0 


30.12 


0.15 


0.98 


6.11 


0.03 


3.5 


3.0 


136.55 


0.36 


9.19 


8.89 


0.07 



Note. — See Table [7] for definitions of quantities. 



TABLE 9 
Encounter Statistics for model Q0.75D1.6 



log(Me/M0) 




A^cnc 


^Poisson 


•^sample 


^, median 


a-r 




[g cm-2] 








[km 


■S-M 




1.5 


0.1 


2.69 


0.04 


0.00 




2.13 


0.01 


1.5 


0.5 


8.01 


0.06 


0.01 




3.23 


0.01 


1.5 


1.0 


12.86 


0.08 


0.02 




3.66 


0.01 


2.0 


0.1 


7.56 


0.07 


0.03 




3.02 


0.01 


2.0 


0.5 


22.05 


0.14 


0.12 




4.46 


0.01 


2.0 


1.0 


34.47 


0.17 


0.18 




4.99 


0.01 


2.5 


0.1 


14.03 


0.10 


0.12 




3.57 


0.01 


2.5 


0.5 


60.88 


0.23 


0.66 




5.25 


0.01 


2.5 


1.0 


132.85 


0.50 


5.88 




6.76 


0.04 


3.0 


0.1 


34.22 


0.19 


1.45 




4.58 


0.04 


3.0 


0.5 


136.45 


0.43 


4.86 




6.53 


0.04 


3.0 


1.0 


275.40 


0.71 


16.94 




8.52 


0.06 


3.5 


0.1 


50.70 


0.16 


2.02 




5.66 


0.04 



Note. — See Table [7] for definitions of quantities. 



TABLE 10 
Encounter Statistics for model GAS 



log{Mc/M0) 


[g cm-2] 


A''giic 


^Poisson 


•^sample 


, .median 

[km s-i] 


Ur 


1.5 


0.1 


2.79 


0.05 


0.01 


2.24 


0.02 


1.5 


0.5 


13.08 


0.12 


0.04 


3.12 


0.01 


1.5 


1.0 


23.42 


0.17 


0.07 


3.78 


0.01 


1.5 


3.0 


53.34 


0.30 


0.20 


4.92 


0.01 


2.0 


0.1 


3.70 


0.05 


0.02 


2.91 


0.01 


2.0 


0.5 


19.62 


0.14 


0.11 


4.29 


0.01 


2.0 


1.0 


35.84 


0.21 


0.23 


5.20 


0.01 


2.0 


3.0 


134.21 


0.51 


1.61 


6.55 


0.01 


2.5 


0.1 


4.97 


0.06 


0.06 


3.59 


0.02 


2.5 


0.5 


25.96 


0.17 


0.46 


5.48 


0.02 


2.5 


1.0 


57.29 


0.27 


0.66 


6.74 


0.01 


2.5 


3.0 


250.26 


0.70 


4.70 


9.27 


0.02 


3.0 


0.1 


5.47 


0.07 


0.13 


4.26 


0.03 


3.0 


0.5 


29.04 


0.18 


0.74 


6.98 


0.03 


3.0 


1.0 


88.33 


0.32 


2.91 


8.53 


0.03 


3.0 


3.0 


318.47 


0.77 


27.18 


10.79 


0.09 


3.5 


0.1 


6.40 


0.06 


0.25 


5.21 


0.04 


3.5 


0.5 


46.90 


0.17 


2.96 


8.54 


0.06 


3.5 


1.0 


99.66 


0.28 


9.34 


9.67 


0.09 


3.5 


3.0 


559.00 


0.79 


16.66 


14.13 


0.03 


4.0 


0.1 


9.62 


0.05 


0.44 


6.25 


0.05 


4.0 


0.5 


63.91 


0.14 


4.20 


9.30 


0.07 


4.0 


1.0 


95.96 


0.22 


3.08 


11.94 


0.03 



Note. — See Table [7] for definitions of quantities. 
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Fig. 2. — Distribution of impact parameters for run Q0.75D2.2, 
with Mc = IQl-^ Mq, Ec = 0.1 g cm-2 (top), and for Mc = lO"'-^ 
Mq, Sc = 1.0 g cni~^ (bottom). Black plus signs show the 1 + 1 
encounters, orange triangles the 1 + 1 + 1 encounters, purple stars 
1 + 2 encounters, and blue squares are the sum of all interactions 
with a well-defined impact parameter. In all cases data points 
show the results of the simulation, with error bars indicating the 
1(T Poisson error, and lines show linear best fits to the data. 



results are similar to those of IParker fc Quan3 (j2012[ ). 
who note that the stripping rat e of planets frorn paren t 
stars decreases with time, and IParker fc Meveil ()2012t ). 
who find that the surface density decreases rapidly after 
one crossing time. 

For a relaxed, bound cluster, the mean number of en- 
counters per star per crossing time (and thus over the 
cluster's entire life for a dispersing cluster) is a function 
of the cluster surface density alone^ and does not depend 
on the mass (jDukes fc KrumhoIS 120121 ). We find that 
this is not the case for unrelaxed fractal clusters. Figure 
[5] shows the number of encounters per solar mass star as 
a function of the cluster mass and surface density; clearly 
the number increases with both mass and surface density. 

We can understand this result by realizing that the 
manner in which our fractal clusters are generated leads 
to an implicit dependence of the surface density on the 
mass of the cluster. This dependence arises because, al- 
though the mean surface density of the cluster averaged 
over its entire face is mass-independent, as the cluster 
mass increases at constant Ec and D the stars become 
packed into smaller and smaller substructures. In the 
Appendix we derive an expression for the effective sur- 
face density Sc,off as a function of Mc and D. Figure |6] 
shows the same data as Figure [SJ but plotted using this 
effective surface density rather than the nominal one. As 
shown in the Figure, the number of encounters is in fact 
nearly independent of Mc and fixed Ec^eff- 

One should regard this result with caution, since it is 
not clear that it remains valid for real clusters, which may 
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Fig. 3. — Distributions of encounter relative velocities in model 
Q0.75D2.2, for the cases (Mc,Sc) = (10l-5,0.1) (top), (103-5,3.0) 
(middle), and (4.5, 1.0) (bottom), where masses are in Mq and Sc 
in g cm^^. Plus signs show data measured from the simulations, 
and lines show the best-fitting Maxwellian distribution. Typical 
reduced x^ values are of order 100, reflecting the poor fits. Notice 
the deviation is largest where the relative frequency of complex 
interactions is higher. 

or may not be truly fractal in their stellar distributions. 
Any attempt to define either Ec or the volume density pc 
for a fractal cluster necessarily requires specifying an av- 
eraging scale over which the quantity is to be measured. 
Ec^cff is best considered as the surface density obtained 
by a process which averages over the initial clumps of 
the substructure, as opposed to the entire cluster. Al- 
ternately, one could envision Ec^ff as being closer to a 
mass- weighted surface density, as opposed to Ec, which 
is an area-weighted surface density. Our result simply 
shows that, if clusters are fractal, then more massive ones 
will produce more encounters than one might guess from 
their surface densities averaged over large scales. This is 
because in a fractal cluster the surface density increases 
as one averages over smaller and smaller scales in the 
vicinity of individual stars. 
It is also interesting to compare our measurements to 
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Fig. 4. — Probability distribution of encounters in time in model 
Q0.75D2.2 for the case Mc = lO^-^ Mq, Sc = 0.1 g cm-^. The 
crossing time for this model is tcr = 0.15 Myr, so the initial peak 
is less than one crossing time in duration. 

the results of a naive analytic estimate. For a uniform, 
spherical, virialized cluster of mass M^ and surface den- 
sity Sc, the expected number of encounters with impact 
parameter b or less in a single crossing time is 



N, 



cnc.cxp 



2Trb^ 



1.2&^i;o 



(11) 



where 63 — b/10^^ AU, Eq — Sc/1 g cm^^, and we have 
used the mean stellar mass from our chosen IMF. Clearly 
the actual number of encounters we measure exceeds this 
value by a large margin. 

3.1.3. Encounter Velocities 

The typical encounter velocity also depends on Mc and 
Tic- Since we found that the distribution of relative ve- 
locities was non-Maxwellian, rather than reporting a ve- 
locity dispersion we instead compute the median v^o of 
the encounters as a function of Mc and Ec. We plot the 
result in Figure [71 The dependence on Mc and Ec is 
somewhat unexpected. For a relaxed cluster, the stellar 
velocity dispersion obeys 



a^ ex (AfcE, 



\0.25 



(12) 



For our non-relaxed clusters, we do see a general increase 
in the relative velocity with mass and surface density, as 
predicted by equation ((T^ . The increase is not as fast as 
expected, however. Figure[S]shows v^a versus Mc at Ec — 
0.1 g cm~^ - effectively a horizontal cut through Figure[7l 
The slope of the best-fit line is approximately 0.15 rather 
than 0.25, and this is true for each of the other surface 
density bins as well. As with the overall non-Maxwellian 
distribution, this slower than expected increase is a result 
of the correlated positions and velocities. 

3.2. Unbound Case (Q1.25D2.2) 

In the case of an unbound cluster, the shapes of the dis- 
tributions of b and Voa are quite similar to those found 
in the base case Q0.75D2.2, so we do not discuss them 
further. The number and median velocity of encoun- 
ters, however, differ from the base case in interesting 
ways. Figure [S] shows the mean number of encounters 
per Solar- type star for the case of an unbound cluster. 
As expected there are fewer encounters in this model 
than there are in the case where at least some of the 
cluster remains bound. Figure [TOl shows the the median 
Woo. Rather surprisingly the velocities tend to actually 
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Fig. 5. — The number of encounters a typical solar mass star can 
expect to experience after 5 crossing times in a fractal cluster, for 
the model Q = 0.75 and D = 2.2 as a function of the cluster mass 
and surface density. 
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Fig. 6. — Same as figure [5] but with the effective surface density 
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Fig. 7. — Median ricx) as a function of Mc and Ec for run 
Q0.75D2.2. 

be slower in this model than in the case of Q = 0.75, de- 
spite the fact that the initial velocities are larger at the 
same Mc and Sc than in the base case. This behavior is 
a result of the initially correlated velocity structure. In 
the base case the cluster tends to relax towards equilib- 
rium, destroying the velocity substructure in the process. 
Once the substructure is destroyed the velocity vectors 
are randomized, producing larger iioo values than during 
the period when the velocity structure is coherent. In 
contrast, an unbound cluster tends to blow apart before 
substructure can be erased, leading us to a lower median 
velocity for those encounters that do occur. 

3.3. High Substructure (Q0.75D1.6) 
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Fig. 8. — Median Vca as a function of the cluster mass at a fixed 
surface density Sc = 0.1 g cm~^ in run Q0.75D2.2. The crosses 
are the simulations results and the line is the best linear fit, which 
has slope 0.15. 
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Fig. 9. — Average number of encounters for a solar mass star in 
an unbound cluster (model Q1.25D2.2) as a function of Mc and Sc 
(top), Scoff (bottom). 
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Fig. 10.— Same as figure[7]but for Q = 1.25 and D = 2.2. 
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Fig. 11. — Same as figure[5]but for D = 1.6 
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Fig. 12.— Same as figure[7]but for D = 1.6. 

The case of high substructure is the most extreme case 
we study in this paper, although it turns out quahtatively 
to be similar to the base case. Figures [TT] and [12] show 
the encounter statistics for this case. Both the number 
of encounters and the relative velocity are higher for this 
case than for the base case. This is consistent with the 
interpretation of D as an increase in the effective surface 
density. However when we plot the number of encounters 
with Eccff the number of encounters in the D = 1.6 
cluster is typically higher than the corresponding point 
in the base case, and the contours are still largely vertical 
rather than horizontal. Thus our model for Ec.off is not 
fully capturing the increase in encounters that occurs 
for very highly substructured clusters containing large 
numbers of stars. 

3.4. Gas Case (Q0.3D2.2) 

Since our Gas runs are extremely subvirial, we ex- 
pect these clusters to relax and destroy substructure ex- 
tremely quickly. Figure [14] shows the temporal distri- 
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Fig. 13. — The distribution of impact parameters for the GAS 
case with Mc = IO'^'^Mq and Sc = 0.5 g cm~^. Here we include 
only events which occur between one and four crossing times (the 
relaxed phase of our model). Black crosses are 1+1 events and 
purple stars are 2+1 events. 

bution of encounters for this model, "which is consistent 
with this expectation. For the first crossing time, there 
is a highly elevated encounter rate as stars fall toward 
the center of the potential well and interact. After this 
they revirialize and the encounter rate drops and be- 
comes roughly constant. Once gas is removed after four 
crossing times, the cluster disperses and the encounter 
rate drops to ve ry small values. 

The results of lParker fc GoodwinI (2012) suggest that 
the distribution function of impact parameters (equation 
([5])) could be significantly altered even in relaxed clusters 
due to the presence of intermediate separation binaries. 
To investigate this possibility, in Figure [13] we shows the 
ensemble distribution of impact parameters for our Gas 
model, considering only encounters that occur at times 
from l—Atc- During this phase the cluster is well-relaxed, 
since it is old enough to have lost most of its initial sub- 
structure, but we have not yet removed the confining gas 
potential. As the Figure shows, the distribution function 
still follows P{b) oc b for 1-1-1 interactions, although some 
evidence of stochastic behavior is observed for the 2+1 
encounters. We have typically seen such deviations for 
2+1 encounters (such as figure [2]), and this is probably 
more due to the difficulty of assigning orbital elements 
when a binary interacts with the single star, as opposed 
to the conditions within the cluster more broadly. 

The number of encounters is shown in Figure 1151 and 
is larger than the number of encounters in the base case 
(shown in Figure [S]), but only slightly - certainly by less 
than the factor of 4 difference in times for which the 
cluster survives before dispersing. Contours of encounter 
number in the (Sc, Mc)-plane are somewhat fiatter for 
the Gas case than the base case, and are much flatter 
in the plane of {T,c,eS,Mc)- Comparing Figures [16] and 
[7] we see that median encounter velocities are larger in 
the Gas case than in the base case. Figure [T7] shows 
the results of taking a horizontal slice to see how the 
median encounter velocity behaves as a function of Mc 
at fixed Sc- We find that the median encounter velocity 
is increasing more quickly with mass in the Gas case 
than the base case, but not quite as quickly as for a fully 
relaxed cluster. 

We can summarize all of these results by saying that 
the Gas case represents a compromise between the case 
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Fig. 14. — The temporal distribution of encounters for a typical 



gas case. This particular model has Mc 



10"^ 



cm ^, so the crossing time is tc 
at 4icr = 0.31 Myr. 



niQ, 



0.5 



: 0.078 Myr. Gas removal occurs 




log M^ [Mq] 




350.00 

315.00 

280.00 

245.00 

210-00 

175-00 

140-00 

105-00 

70-00 

35-00 

0.00 



2 3 

log M^ [Mq] 



Fig. 15. — Same as figure[5]but for the gas case. 

of a fully relaxed cluster and the substructured clus- 
ters we have considered thus far; the first crossing time, 
which produces a significant fraction of the encounters 
due to the elevated encounter rate during the relaxation 
phase, looks much like the substructured clusters. After 
a crossing time the stars revirialize, and the evolution 
from that point until gas expulsion looks like that in a 
dynamically relaxed cluster. Because a significant frac- 
tion of the encounters occur during the early, subvirial 
relaxation phase, the overall number of encounters and 
encounter velocity distribution ends up being intermedi- 
ate between the substructured and fully-relaxed cases. A 
lifetime of four crossing times before gas expulsion is not 
long enough for the overall statistics to be dominated by 
the relaxed phase rather than the unrelaxed one. 

4. DISCUSSION AND IMPLICATIONS 
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Fig. 16. — Same as figure[7]but for the gas case. 
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Fig. 17. — The log of the median Vao versus the log of Mc at 
Sc = 0.1 g cm"-^. The red crosses are our data points with the 
best-fit line drawn in black. The slope in this case is 0.17. 

4.1. Error Analysis 

How certain are the values of A^cnc derived above, given 
the number of simulations we have run and the number 
of encounters they produce? This question requires some 
care. There are two distinct sources of error, and each 
dominates in a different regime of our simulations. One 
source of error is simply counting statistics on the total 
number of events at a given (AfcSc)- At low values of 
Mc and Sc, even when we have a large number of runs, 
the total number of events recorded over all simulations 
may be small, producing a significant statistical error. 
The second source of error arises from our limited sam- 
pling of all possible realizations of fractal clusters at a 
given {Mc, Sc)- At large {Mc, Sc) the number of events in 
a given run may be quite large, producing small Poisson 
errors, but the numbers of events may be quite different 
for different realizations at the same {Mc, Ec). For exam- 
ple, we might have three realizations that produce 8,000, 
10,000, and 12,000 events, respectively. In this case the 
Poisson error on each of these numbers is of order 100, 
much smaller than the difference between them, indicat- 
ing that our error is dominated by our limited sampling 
of possible clusters at a given {Mc,T,c)- A reasonable 
estimate of the error in this regime is the standard de- 
viations of the mean values for each run, neglecting the 
Poisson errors on each individual run. To interpolate 
between these two limits, we take the total error to be 
the quadrature sum of these two types of error. This is 
approximate, but should be roughly correct. 

To make the above analysis precise, let A^run be the 
number of simulations run for a given set of parameters 



{Mc,I^c,Q,D). Let Si and Ei be the number of Solar 
mass stars and encounters in the i"' realization, respec- 
tively. Then we define 



Wru 



r 
events " 


= >>. 




j=i 




iV,u„ 


N,-- 


= >>. 



so that 



iV™n = 



i=l 



iV.- 



TV* 



(13) 
(14) 

(15) 



is our best estimate of the number of encounters per Sun- 
like star. From these definitions we can express the error 
on Ncnc as 



rAr2 _ "-''events 2 



7V2 



Poisson ' sample 



cr„ 



where 



fPo 



vn:. 



TV. 



(16) 



(17) 



is the error due to counting statisticeQ and 



f^samplc 



r::T{E^-Ey/N.. 



1/2 



iV. 



(18) 



is the error introduced by a lack of ability to fully sam- 
ple the parameter space by having enough realizations of 
initial clusters. Here E is the mean number of events per 
run. We report crpoisson and CTsampio separately in Tables 
[7] - [ini along with the total relative error 



SN, 



N,. 



4.2. The Sun's Birth Environment 



(19) 



One of the primary applications of the statistics we 
have measured is constraining the environment in which 
the Sun was born. To do so, we make use the the 
velocity-dependent cross sections fo r perturbing the or- 
bits of the outer planets measured bv lDukes fc Krumholzl 
(.201 ^ from their s imulations. In deriving these values 
[bukes fc KrumhoI3 implicitly integrated over the distri- 
bution of impact parameters under the assumption that 
this distribution follows p{b) oc b, and we have found that 
this is generally a good assumption. 

If <^i{voo) is the cross section for a particular event 
to occur for stars approaching with a particular relative 
velocity Wqo , the expected number of occurrences of that 
event in a cluster within which there are expected to be 
A'onc encounters with impact parameter less than 6niax = 
1000 AU is 



A„; 



^enc J^ 0-i{Vao)VooP{Voo) dVo 



max 



J^ VooP{Voo) dVa 



(20) 



■^ Note that this expression assumes the Gaussian limit, and 
therefore becomes invalid when A^ovcnts J; 10. 
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Assuming the events are Poisson in nature (i.e. that they 
are independent), the probabiUty of an occurrence is sim- 
ply 

P, = l-exp(-A,0. (21) 

Following iDukes fc KrumhoI3 (|2012f ). we examine the 
possibihty that a close encounter with a passing star 
might excite one of the Jovian planets to a highly ec- 
centric (e > 0.1) orbifl To eva luate Aj,, we combin e the 
measured values of ai{voo) from lDukes fc Krumhold with 
velocity distributions p{voo) obtained in the previous sec- 
tion. 

Figure [18] shows the probability of exciting a Jovian 
planet to eccentricity e > 0.1 in a cluster with Q — 0.75, 
D = 2.2, as a function of Mc and Sc- Overall our con- 



clusions are consistent with those of IDukes fc Krumhold 
()2012|) : even for clusters of quite high masses and sur- 
face densities (even up to Ec,eff = 20 g cm~^), it is ex- 
tremely unlikely that a Sun-like star would have an en- 
counter with another star close enough to significantly 
perturb the orbit of a planet like Neptune. However, 
we also find that the probability of excitation is inde- 
pendent of Mc at fixed Sc^eff, but that it increases with 
M r at fixed E^. This i s the opposite of the trend found 
bv IDukes fc Krumhold (2012), and the difference is easy 
to understand: for a relaxed cluster, as assumed by 
IDukes fc Krumhold . the number of encounters is inde- 
pendent of Mc at fixed Ec, while the typical encounter 
velocity increases with mass as Mc'"^^, reducing the effec- 
tive cross-section per encounter. For unrelaxed clusters, 
on the other hand, we find that the number of encoun- 
ters increases with Mc and fixed Ec, and that the typical 
encounter velocity increases with Mc more slowly than 
is the case for a relaxed syst em. These two chan ges re- 
verse the trend predicted bv IDukes fc Krumhold . How- 
ever, while the direction of trend with mass is reversed, 
the trend remains rather weak. As a result, the quali- 
tative conclusion that encounters in dispersing clusters 
should have no significant impact on the Solar system, 
even if those clusters are quite massive, continues to hold. 

In the unbound case we find similar results to the base 
case. The combined results of a significant drop in the 
number of encounters, with a slight decrease in the me- 
dian relative velocity leads to a slightly lower excitation 
probability. Figure [12] shows the probability of exciting 
Neptune's orbit. Qualitatively the results are similar to 
those in figure [18] The same is true for the high substruc- 
ture case, model Q0.75D1.6, shown in Figure [20] Com- 
pared to the base case, the effect of increased number of 
encounters is much stronger than the increase in relative 
velocities, leading us to higher excitation probabilities. 
Again the shape is similar, but slightly more extreme 
than the base case. The overall probability remains low 
at low masses, but seems likely to become significant at 
masses above ~ 10^ Mq, provided that nominal (as op- 
posed to effective) surface densities remain roughly con- 
stant. Whether such highly substructured, massive clus- 
ters exist in nature is uncertain. 

The most interesting effects are seen in the Gas case. 

^ Alternately, since the Jovian planets are likely not fully formed 
during the early phases of evolution we are considering, we can 
regard these probabilities as the describing the chances that an 
encounter with another star might severely perturb the protoplan- 
etary disk at the radii where the Jovian planets are located now. 





Fig. 18. — The log of the probability that a Jovian planet is 
excited to an eccentricity e > 0.1 as a function of Sc (top 4) and 
Sc.cff (bottom 4) and Mc, for model Q = 0.75, D = 2.2. Within 
each set of four we have Jupiter (top left), Saturn (top right), 
Uranus (bottom left) and Neptune (bottom right). 

Due to the destruction of substructure because of the 
subvi rial nature of the cluster, this case approaches that 
of D ukes fc KrumhoI3 ()2012D . Figure [21] shows the prob- 
ability of exciting Neptune's orbit in such a cluster. The 
increase in probability of disruptive events with cluster 
mass is very weak in this case, and when we consider 
the effective surface density we are able to reproduce 
the trend of IDukes fc Krumhold . Namely, at constant 
(effective) surface density, we find that the probability 
of a disruptive event actually decreases with the clus- 
ter mass. This implies that there is no dynamical limit 
on the number of stars in the stellar birth cluster. The 
overall disruption probabilities are slightly higher than 
those of the base case, but this is to be expected since 
the cluster undergoes an initial period of highly elevated 
encounter rate before relaxing and then dispersing. 

We can also estimate the errors on these probabilities 
as follows. The error on the probability of any event 
(e.g. excitation of a Jovian planet) is given by 



SPi = exp{~A,)SAi. 
With a little algebra one may show that 



i6A,r = A^ 



SNc. 



Nc: 



h 



(Sh 

\l2 



(22) 



, (23) 



where /i and I2 are the integrals in the numerator and 
denominator of equation (j20p . respectively. 

The errors on /i and I2 may be obtained by differen- 
tiation of the Riemann sum approximating the integral: 



iShf 



E2 2 2 



(Sa,)^ 
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Pt 
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(24) 
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Fig. 19. — The log of the probability of exciting Neptune to an 
eccentricity e > 0.1 as a function of Sc (top) and Sc,cff (bottom) 
and Mc, for model Q = 1.25, D = 2.2. Notice that the probabilities 
are lower than those in the unbound case. The other Jovian planets 
have similar plots to these. 
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Fig. 20. — Same as figure [Tol but in the case of high substructure 
(D = 1.6, Q = 0.75). The other Jovian planets have similar plots, 
but with a lower excitation probability. 



Fig. 21. — Same as figure FlQl but for the case where we include 
an external (gas) potential for 4 crossing times, and then allow the 
cluster to disperse. 

For the sake of compactness of the expression we omit 
the oo subscript on the velocity. I2 is given by a similar 
expression, except that the cross section is absent. Given 
the errors on A^cnc (see Section [4T]) . /i, and I2, we can 
compute the overall errors on our probabilities from (|22p . 
To keep this analysis as simple as possible we assume 
that the relative errors on the integrals are very small, 
so that the relative error on the number of encounters 
dominates. Then we have 



SAi = A. 



SN,-_ 



N,. 



(25) 



Typically A^ <C 1 (all cases except one have ANcptunc < 
10^^, and all of the other Jovian planets have smaller 
values) and in this case we can expand equation (j22p to 
first order in Ai to obtain 



dFi « Ai— = AiCTr 



(26) 



Since typical values of A^ are 0(10 ^) or smaller, and 
typical values or ar are also O(10~^) or smaller, this 
means that the absolute error on the percentages calcu- 
lated through this method are less than 0.1%, and the 
relative error, of order ar, is at most ^ 10%. Thus we 
can be fairly confident in our major results. 

Finally, there a few additional sources of error that 
we mention here, but do not quantify explicitly. The 
IMF used in lDukes fc Krumholj (j2012t ) is slightly differ- 
ent than ours, and yields a mean stellar mass of approx- 
imately O.2M0 instead of O.SOMq. This implies that 
we would expect the value of the probabilities to be 
slightly higher, although the shape of the constant prob- 
ability curves would not change (figure ITS)) . Similarly, 
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the fact that there is shght deviation from the distri- 
bution P{b) ex b, especially at high Ec, is a source of 
error. At low surface density, however, this error should 
be minimal. To obtain an entirely correct result for Aj 
one should repeat the calculations of Duk es fc Krumhold 
(|2012 ') using the same IMF and distribution of impact 
parameters as found in these simulations. 

4.3. Free-Floating Planets 

Ever since their discovery bv lSumi et al.l (|2011l ). there 
has been considerable debate about the origin of plan- 
etary mass objects that are either completely unbound 
or very distant from their parent stars. Models for the 
origin of these objects include planet-planet scattering 
leading to ejection in a young planetary system (e.g. 
iNagasawa fc Idal[2011l : lBolev et a l. 2012), escape caused 
by mass loss dur ing late stages of stellar evolution (e.g. 
iVeras et al.ll2011[). direct form ation from the interstellar 
medium (jStrigari et al.l 120121) . and d ynamical stripping 
of planets from stars in . clusters (e.g. iBolev et al.l 120121 : 
IVeras fc RavmondllMl fOvelar et al.ll2012D . While a full 
discussion of this topic is beyond the scope of this paper, 
we do note that our simulations provide new insight into 
the latter mechanism. 

Our results from the previous section show that signif- 
icant orbital disturbance for a planet at 30 AU, the dis- 
tance of Neptune, is likely to be a relatively rare event. 
Only for clusters whose masses exceed ~ 10^ Mq and 
are very highly substructured {D = 1.6) do orbital exci- 
tations become common. Complete ejection will be even 
rarer. This suggests that, unless the planet formation 
process is capable of producing large numbers of plan- 
etary mass objects at radii significantly beyond 30 AU 
by internal mechanisms (e.g. planet-planet scattering), 
cluster stripping cannot be a significant source of free- 
floating planets. A more precise and quantitative ver- 
sion of this statement may be obtained by combining the 
orbital element distributions we have obtained here with 
population synthesis models for star clusters and planets. 

Finally, it is interestin g to compare the results of our 
simulations with those of ' Parker fc Ouana (|2012( ). They 
find large orbital excitations only in their Q = 0.3 case, 
which includes no gas potential term. This will lead to a 
tightly bound final cluster. This fact together with their 
long simulation times (10 Myr, or ~ 30 crossing times) 
means that it is not surprising that there are a large 
fraction of planets becoming unbound. In addition, they 
consider orbital excitations of planets from stars of any 
mass, and by number the majority of stars are smaller 
than the 0.8 — I.2M0 mass range that we consider. The 
binding energy of planets in such systems will be smaller 
than is typical of the planetary systems we consider. 

The closest match in parameters be tween our sim- 
ulations and those of iParker fc Quanzj is between our 
Q = 0.75,1? = 2.2 and their Q = 0.7, D = 2.0 runs. 
Their figure 8 corresponds roughly to clusters with Mc = 
IO^-IO^Mq and E^ = 0.1 g cm'^. They find that - 10% 
of planets at 30 AU are stripped from their parent star 
in this case, whereas we find typically that ~ 2 — 3% 
of orbits are excited in such a case. These differences 
are most likely due to the issue of binding energies m en- 
tioned above. For the IMF used by IParker fc OuanH and 
their assumption that planets are equally likely to occur 
around stars of any mass, only ~ 20% of planets or- 



bit stars of mass > 1 Mq, while ~ 40% are born around 
stars with mass < 0.5 Mq. Thus, it is not surprising that 
they find a higher rate of orbital perturbation and strip- 
ping. Which model will more accurately predict numbers 
of free-floating planets is unclear, due to the lack of reli- 
able estimates of planet fractions around stars with mass 
< IMq. 

5. SUMMARY AND CONCLUSIONS 

We have conducted a series of simulations of the evo- 
lution of dispersing young star clusters that possess a 
high degree of substructure. These simulations cover a 
wide range of masses (30 — 30,000 Mq), surface densities 
(Ec = 0.1 — 3.0 g cm~^), dynamical states (supervirial 
but bound, unbound, and subvirial but then unbound 
by gas expulsion), and degrees of substructure (fractal 
dimension D = 1.6 and 2.2). These parameters are cho- 
sen to reproduce the range of properties for young star 
clusters suggested by both observations and star forma- 
tion simulations. We provide tabulated distributions of 
number of encounters, impact parameters, and relative 
velocities as a function of these properties. These may be 
used as inputs in population synthesis models of planet 
formation. 

Our calculations produce a number of interesting con- 
clusions. First, during the '^ 1 crossing time it takes 
for the initial substructure to be erased by dynamical 
interactions, the number of encounters is significantly el- 
evated compared to what one would expect for a relaxed 
system. The amount of elevation depends the amount 
of substructure and other cluster properties. Regardless 
of its strength, though, the enhancement is transient. 
Either the substructure dissolves (if the cluster is not 
confined or mildly bound), or the cluster disperses (if 
it is strongly unbound). Thus for moderate degrees of 
substructure the overall enhancement in the number of 
encounters is only a factor of a few for clusters that do 
remain bound for some period before gas dispersal. Only 
if the gas has extremely strong substructure (fractal di- 
mension D = 1.6) is the enhancement larger. 

Second, early in the evolution of a substructured clus- 
ter, before the cluster has dispersed or the substructure 
has been erased, the distribution of encounter impact pa- 
rameters is not far off from the expectation for a relaxed 
cluster, but the distribution of velocities is significantly 
non-Maxwellian. Because this early phase contributes an 
appreciable fraction of all encounters even in clusters that 
remain bound for four crossing times, the overall distri- 
bution of encounter velocities is non-Maxwellian even in 
such clusters. Compared to a Maxwellian, our clusters 
show both a sharper peak at moderate encounter veloci- 
ties, and a longer tail extending to higher velocities. The 
overall median velocity increases with cluster mass more 
slowly than one would expect in a relaxed cluster, and 
scales with the virial velocity ratio. Clusters with larger 
virial ratios, and thus larger velocity dispersions, actu- 
ally tend to produce lower median encounter velocities 
because they are less effective at dissolving the velocity 
substructure and randomizing stellar relative velocities. 

Third, even with the enhanced encounter rates that we 
find, we conclude that planetary systems or protoplane- 
tary disks around stars in dissolving clusters are unlikely 
to experience significant dynamical perturbations from 
other stars in the cluster, at least for planets that are 
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TABLE 11 

Number of Generations Versus Cluster Mass 



D 


1.5 


2.0 


log(Mc/M0) 
2.5 3.0 3.5 


4.0 


4.5 


2.2 
1.6 


3.02 
3.99 


4.00 
4.98 


4.99 5.42 6.01 
5.99 7.01 8.00 


7.00 


7.98 



Note. — The average number of generations re- 
quired to generate a cluster of mass Mc, for D = 2.2 
and D = 1.6. 

within tens of AU of their parent star. Such planets are We thank D. Dukes for help running simulations and 
simply too tightly bound and have cross sections that analyzing the data, and S. Aarseth for assistance with 
are too small for many of them to be disturbed. This re- NB0DY6. MRK acknowledges support from the Alfred 
mains true even in our most highly substructured cases, P. Sloan Foundation, the NSF through grant CAREER- 
which produce the largest number of encounters, up to 0955300, and NASA through Astrophysics Theory and 
cluster masses of 10^'^ M©. This means that there is Fundamental Physics Grant NNX09AK31G, and a Chan- 
no dynamical constraint on the size of the Sun's parent dra Space Telescope Grant, 
cluster, and that cluster stripping is unlikely to be an 
important contributor to the population of free-floating 
planets in the Milky Way. 

APPENDIX 
A. DERIVATION OF AN EFFECTIVE SURFACE DENSITY FOR FRACTAL STAR CLUSTERS 

We can define an effective surface density Sc.off by considering the process by which the fractal cluster is generated. 
Let the cube out of which the fractal is built have a volume given by Vq, with a characteristic radius r^. The process 
of building the fractal can be thought of as removing chunks of the volume from this initial value. Since the first 
generation particles are always parents, the volume loss starts at the second generation. The volume lost after the 
second generation is, on average, Moss, 2 = Vo(l — 2^^'^). The volume remaining after the second generation is then 
K-omain,2 = Vq - Moss,2 = Vq2°'''^ . Similarly, the volume lost after the third generation is Moss,3 = Vrcmain,2(l - 2^"^), 
the remaining volume is V^emain,3 — Vb(2^~^)^, and so forth. The volume remaining after g generations is simply 



''romain.g — ''O^ ■ \-^^) 



We can thus define an effective radius by 



r,,,g=r,2i(^-3)(,-i)^ (A2) 

Here we have implicitly assumed that we have taken enough generations of the fractal that when we make the distri- 
bution of positions spherical, that the volume loss remains the same. By using dimensional scaling arguments we see 
that the effective surface density is then 

I]e,cff(/?,5)-S22-i(^-3)(9-i). (A3) 

where YP^ is the surface density of a cluster with constant density of stars and radius Vc- The effective initial surface 
density therefore depends on how we choose the number of generations for our fractal, which in turn is determined by 
the mass of the cluster, since the number of generations required is determined by the condition that there be enough 
potential sites to accommodate the number of stars in the cluster. We approximate this condition by 

.W-^ + l + =s.p), (A4) 

where N = Mc/fh is the number of stars in the cluster and S2{D) is 1 if I? < 2 and otherwise. This choice comes 
from the fact that at _D = 2, Pparcnt = 0.5. The average number of generations required to generate the fractal as 
a function of the mass is displayed in table [TTJ We obtained this result by averaging over a random sample of 100 
different fractals for each mass bin. 
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